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SECTION  1 
INTRODUCTION 

Sandwich  composites  constitute  an  important  class  of 
materials  in  the  aerospace  and  building  industries  due  to  their 
potential  to  provide  high  resistance  to  loading  at  a  relatively 
low  weight  penalty.  Typical  uses  of  sandwich  materials 
presently  include  wing  skins  and  control  surfaces,  fairings, 
shelving,  cargo  doors,  helicopter  blades,  and  prefabricated 
panels  for  building  construction. 

As  the  application  of  these  advanced  materials  increases, 
sophisticated  analysis  techniques  become  increasingly  important 
for  use  in  the  formulation  and  qualification  of  practical 
designs.  The  necessary  analytical  methods  must  be  readily 
accessible  to  the  designer,  while  remaining  quite  general  in 
scope.  To  achieve  the  objective  of  generality,  computational 
methods  must  be  employed,  the  most  flexible  of  these  being 
techniques  based  upon  finite  element1  discretization. 

The  development  of  finite  element  analysis  techniques 

applicable  to  sandwich  constructions  is  complicated  by  the 

nature  of  the  response  of  most  sandwich  layups.  The  relatively 

thin,  high-modulus  face  sheet  layers  resist  loading  through 

inplane  stresses,  and  thus  are  represented  most  effectively 

with  plate  or  shell  type  elements.  Sandwich  core,  which  is 

typically  much  thicker  and  more  flexible  than  the  face  sheets, 

deforms  principally  in  transverse  shear  modes,  and  therefore, 

is  modeled  most  appropriately  by  three-dimensional,  shear- 

deformable  elements.  The  incompatibility  of  these  two  classes 

of  finite  elements  (shells  and  solids)  is  well-known,  and  a 

departure  from  the  traditional  elements  of  these  types  is 

clearly  necessary.  A  number  of  specialized  formulations  have 

been  proposed  to  deal  with  the  modeling  problems  of  sandwich 
2—6 

materials  .  However,  none  of  these  has  been  developed  in 
sufficient  generality  to  provide  truly  comprehensive  modeling 
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capabilities  for  sandwich  panels  having  arbitrary  curvature, 
multicore  construction,  transitions  to  other  types  of  structure, 
edge  closeouts,  and  other  troublesome  geometrical  features. 
Applications  to  nonlinear  sandwich  analysis,  including  finite 
displacements,  large  rotations,  and  material  nonlinearities 
are  typically  further  restricted  (flat  panels,  single-core 
sandwich,  or  rectangular  shapes)  due  to  the  complexity  of 
the  formulation  and  long  solution  times. 

The  present  report  documents  a  finite  element  approach 
for  the  analysis  of  sandwich  structures  of  arbitrary  geometry, 
which  may  include  stiffening  members,  closeouts,  and  connections 
to  other  three-dimensional  structural  components.  Static 
response  to  mechanical  and  thermal  loading  may  be  computed  for 
either  linear  or  nonlinear  deformations.  In  the  case  of 
nonlinear  analysis,  arbitrarily  large  displacements  and 
rotations  are  permitted,  in  addition  to  material  nonlinearities 
(plasticity) .  Linear  natural  frequency  and  normal  mode 
calculations  can  also  be  performed. 
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oC.CT  I  ON  2 

THEORETICAL  DEVELOPMENT 


The  basis  of  the  present  analysis  of  sandwich  materials 
is  the  incremental  equilibrium  relation  for  a  general,  three- 
dimensional  continuum  experiencing  large  displacements  and 
materially  nonlinear  response.  Each  of  the  finite  elements 
(face  sheet,  core,  stiffener,  solid)  described  in  Section  3 
is  derived  from  this  general  set  of  governing  equations  by 
first  making  appropriate  specializations  and  then  applying  the 
procedure  of  finite  element  discretizations. 

For  finite  element  application,  the  appropriate  form  of 
the  equilibrium  equations  is  expressed  in  terms  of  the  principle 
of  virtual  work.  In  this  Section,  the  incremental  form  of  the 
principal  of  virtual  work  is  obtained  for  a  general,  three- 
dimensional  continuum.  The  Lagrangian  description  of  motion 
is  used  throughout:  that  is,  all  displacements,  strains  and 
stresses  are  expressed  in  terms  of  the  original  (undeformed) 
configuration  of  the  body.  The  appropriate  constitutive 
relations  for  elastic  and  elastic-plastic  material  behavior 
are  also  developed. 

2.1  INCREMENTAL  EQUATIONS  OF  MOTION 

For  the  purpose  of  obtaining  an  incremental  description 
of  motion,  three  configurations  of  a  general  structure  are 
considered : 

Configuration  C®  :  the  initial  state  of  the  body, 

used  as  a  reference  state 

Configuration  C1  :  an  intermediate  state,  at  which 

a  solution  for  the  structural 
response  is  assumed  tc  be  known 

2 

Configuration  C  :  a  subsequent  configuration, 

removed  from  state  C^  by  a  single 
increment  of  loading. 
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These  three  states  of  deformation  are  shown  in  Figure  1.  In 

an  incremental  nonlinear  analysis,  the  objective  at  each  step 

2 

of  the  solution  is  to  compute  a  configuration  such  as  C  ,  given 
the  previous  solution  (state  C^)  and  the  incremental  values  of 
mechanical  and  thermal  loading. 

In  the  initial  state  C®,  the  position  of  an  arbitrary 

point  P  (Figure  1)  is  denoted  by  X.;  i=l,2,3.  In  reaching 
1  1 

configuration  C  ,  the  point  P  moves  to  a  new  position  whose 
coordinates  are  x^,  where 


xi  -  Xi  +  lui 


and  .u.  represent  the  displacements  in  state  C  .  The  state  of 

ii  i 

strain  at  state  C  is  measured  by  the  Green-St.  Venant  strain 

^  7 

tensor 


lei j  “  2(lUi,j  +  lUj,i  +  luk,  i  luk, '  (2) 

Observe  that,  in  Equation  2,  differentiation  is  performed  with 
respect  to  the  initial  coordinates  X^;  that  is, 

lui,j  "  3lui/3Xj*  (3) 

The  appropriate  measure  of  the  stress  corresponding  to 
Equation  2  are 

,s,  .  =  /G  ,n.  .  (4) 

1  lj  1  lj 

.  7 

where  .11.  .  are  the  symmetric  Piola-Kirchof f  stresses  in 
1  1  _ 

configuration  C  ,  and  /G  is  the  determinant  of  the  Green's 
deformation  tensor 
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Similarly  in  configuration  C  ,  denote  the  position  of  the  point 
P  by  coordinates  z^;  i*l,2,3.  Then 

z^  ■  +  Au^  (6) 

in  which  Au.  are  incremental  displacements  separating  the  states 
1  2  1 

C  and  C  .  The  strain  and  stress  are 

2*i  j  '  l*ij  +  “ij  ,7) 

and 


2*1 j  '  1*11  +  4\j 


respectively. 


Next  consider  the  principal  of  virtual  work  in 
1  2 

configurations  C  and  C  ,  which  are  written 


\  i9ijxk,iSuk,jdv  -  f  i?iSoi" 


-s.  .  z.  .  6u.  .dV  =  «t.6u.dA 

J  2  13  k,i  k,  j  J  2  l  i 


Here  V  represents  the  original  volume  of  the  body  (in  state 
0  ° 

C  )  ,  3V  denotes  the  corresponding  boundary  surface,  and 

-  1  2 
lfci'  2fci  are  Prescrike<*  surface  forces  in  states  C  and  C 

respectively.  It  is  important  to  notice  that  the  virtual 

displacements  6u^  may  be  any  kinematically  admissible 

displacements,  and  that  they  are  measured  in  each  case  with 

respect  to  the  initial  configuration.  Thus  by  subtracting 

Equations  9  and  10  it  is  possible  to  write 


(,s,  ,r.  .  -  ,g.  ,x.  ,  )Au. 


2°ij‘ic,i  -  i“irk,i'-k,r*  -  i  i 


L* 


6u .  dA 


(11) 


where  St^  are  the  increments  in  prescribed  surface  traction, 

Ki  “  2*i  “  lV  (12) 

Using  the  incremental  stresses  defined  in  Equation  8,  Equation 
11  can  be  arranged  in  the  form 

f  <i*ij4uk.i  ♦  4,ijzk,i)  “uk,jdv  '  j  IV4uidA  <13) 

V  _  3  V 


Note  that  the  arbitrary  virtual  displacements  6u^  have  been 
replaced  by  incremental  virtual  displacements  6Au^,  which  is 
clearly  admissible  provided  both  sets  of  virtual  displacements 
satisfy  the  imposed  kinematic  boundary  conditions  of  the  problem. 

Equation  13  expresses  the  equality  of  internal  and 

external  virtual  work  performed  during  the  increment  of 

1  2 

deformation  between  configurations  C  and  C  .  Using  Equation 
6,  the  equation  obtained  above  can  be  expressed  more 
conveniently  in  terms  of  the  incremental  strains  Aeij  (Equation  7) . 
Since  the  incremental  stress  tensor  As^  is  symmetric. 


As  .  ,  z.i.  <SAu.  . 
i j  k  l  k»  3 


=  2  4sij  (Zk,i44uk.j  +  1141 


and  thus 


4*ijzk,i51uk,j  ■  4sijS4«ij-  (15> 

The  incremental  principle  of  virtual  work.  Equation  13,  then 
becomes 


7 


( 


f  (Aalj<A*1J  ♦  1«ljA»k,i«A“x,j)aV  * 
ov 

|  SE^AAUjdA. 
o*v 

For  later  use,  it  is  useful  to  separate  those  terms  of  the 
incremental  strain  Ae^  which  are  linear  and  nonlinear  in  the 
incremental  displacements  Au^. 

Define 


4eij  *  4Eu  +  4"i3 


a 


in  which 

4Eij  ’  \  (4ui,j  +  4uJ.i  *  iuk,i4uk,:)  +  1^, jAuk, i>  (] 

4nij  "  5  luk,i4“k.j '  ( 

Equation  (17)  becomes 

|  (Asi j6Aei j  +  jS^iAn^dv  » 
oV 

|  At^SAu^dA  ( 

o ™ 

The  material  stress-strain  relation  used  here  assumes 
that  the  incremental  stresses  and  strains  can  be  related  in  a 
linear  fashion. 


Asij  "  Dijkl(Aeij  "  akiAT) ’ 


(2 


8 


Here  are  the  thermal  expansion  coefficients  and  AT  is  an 
incremental  change  in  temperature  at  a  point.  Combining 
Equations  20  and  21  yields  the  final  form  of  the  incremental 
principle  of  virtual  work, 

f  (Dijkt(4eia  -  +  i'ijMiijjav  - 

V 

o 

|  Ati6AuidA.  (22) 

o3V 

Equation  22  is  fully  three-dimensional,  and  specific  forms 
which  are  used  to  derive  the  individual  finite  element  types 
are  discussed  in  Section  3.  The  iterative  solution  of 
Equation  22,  which  is  nonlinear,  is  outlined  in  Section  4. 

2.2  CONSTITUTIVE  RELATIONS  (ELASTIC  MATERIALS) 

For  elastic  materials,  the  incremental  stress-strain  law 
of  Equation  21  is  precisely  the  same  as  the  relation  between 
total  stresses  and  strains, 

sij  -  Dijkt<V  -  «k*T,<  (23) 

where  T  is  interpreted  as  the  temperature  change  from  a  fixed 
reference  temperature.  In  matrix  form.  Equation  23  can  be 
written  as 


s  -  D  (e  -  oT) 


(24) 


in  which 

®  "  L®11  °22  s33  S23  S13  S12-f 
e  -  Len  e22  e33  e23  e13  e12-J 


(25) 

(26) 
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and  D  is  a  coefficient  matrix  of  order  six. 

In  the  case  of  an  isotropic  material,  D  is  completely 
determined  by  the  values  of  Young's  modulus  E  and  the  Poisson's 
ratio  v, 

-v)  v  v 

v  (1-v)  v 

v  v  (1-v) 

0  0  0  ( 

0  0  0 

0  0  0 

and  a  depends  only  upon  the  scalar  coefficient  of  thermal 
expansion ,  ct : 

«T  “  La  a  a  0  0  °J  •  (29) 

For  orthotropic  materials,  the  elastic  properties  are 
direction-dependent,  and  a  stress-strain  law  must  first  be 
formulated  with  respect  to  the  preferential  axes  of  the  material 
and  then  transformed  to  global  coordinates.  In  the  material 
coordinate  system  (denoted  by  a  subscript  m) ,  Equation  25 
becomes 


D  (e  -  a  T) . 

-itl  -Iu  -In 


The  matrix  Dm  is  a  function  of  nine  independent  material 
constants:  the  extensional  moduli  E^E^Ej,*  the  Poisson's 
Ratios  v23,  v13,  v12;  and  the  shear  moduli  G23,  G13,  G^2. 


The 


1 


nonzero  upper  triangular  elements  of  D^,  which  is  symmetric, 


are  given 

by8 

dll 

=  (1  - 

',23V32>VN 

d12 

=  (V12 

+  V32V13)VN 

d13 

=  (V13 

+  V12V23,E3/N 

d22 

=  (1  - 

v13V32)VN 

d23 

"  (V23 

+  ''21u13)VK 

d33 

=  (1  - 

v12''21>VN 

d44 

=  G23 

d55 

=  G13 

d66 

=  G12 

(31) 


in  which 


N  ’  1  '  y12V21  '  y23V32  -  V31V13  '  2u21V32V13' 


(32) 


It  should  be  noted  that  the  material  properties  in  Equations 
31  must  satisfy  certain  constraints  in  order  to  be  valid 
from  symmetry  and  energy  considerations.  These  constraints 

g 

can  be  summarized  as  follows  : 


'L 


>  0 

(33) 

>  0 

(34) 
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(i-vijV  > 0 

(no  sum) 

(35) 

!ii  - 

Ei  Ej 

(36) 

N  >  0. 

(37) 

With  respect  to  the  global  coordinates  of  a  body,  the 
stress-strain  relation  in  Equation  30  becomes 


S  =  T  D  T*  (e  -  T^aT). 

~  ~s~m~e  -  -e  ~m 

(38) 

Here  T  ,  T 
-  s  ~e 

between  the 

define  the  transformations 

material  and  global  axes: 

for  stress  and  strain 

s  -  Ta 

~  ~  s  ~m 

(39) 

e  *  T-  • 

~  e  -m 

(40) 

A  comparison  of  Equations  23  and  38  yields  the  correct 

transformations  needed  for  D  and  a  ,  namely 

-m  -m 

D  -  T  D  T*  (41) 

-  ~  s  ~m  ~e 

a  »  T-^  a  .  (42) 

~  ~  e  ~m 


2.3  CONSTITUTIVE  RELATIONS  (ELASTIC-PLASTIC  MATERIALS) 

The  elastic-plastic  material  law  considered  in  the 

present  development,  which  is  applicable  to  initially  isotropic 

,  9 

materials,  is  based  upon  the  von  Mises  yield  criterion  , 


F  (s 


-  g  -  _ 

i  jSi  j 


(43) 
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in  which  sT..  represent  the  deviatoric  stresses 


8ij  “  8ij  "  38kk5ij 


2 

and  k  represents  the  diameter  of  the  yield  surface.  For  this 
yield  function,  the  associated  plastic  flow  rule  is 


or,  from  Equation  43, 


(44) 


(45) 


That  is,  the  individual  components  of  the  plastic  strain  rate 
are  proportional  to  the  corresponding  deviatoric  stresses. 

The  consistency  condition  is  used  in  the  form  noted  by 
Hunsaker10  to  obtain  a  computationally  effective  procedure. 
Thus 


3F 


3s 


ij 


Sil 


H 


IP 

3®ij  ij 


(46) 


The  parameter  H  characterizes  the  strain-hardening  slope  of  the 
stress-strain  curve,  and  in  the  one  dimensional  case  reduces  to 


H 


E  E„ 


E-E„ 


(47) 


in  which  ET  is  the  tangent  to  the  uniaxial  stress-strain  curve. 

Finally,  assuming  an  additive  decomposition  of  the  strain 
rates  into  elastic  and  plastic  components  gives 


'ij  “  Eijk£(ekl  “  ek£K 


(48) 
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Here  E^..^  is  the  modulus  tensor  of  the  material  in  its  elastic 
state.  Combining  Equations  44,  46  and  48  permits  a  solution 
for  X  in  terms  of  the  total  strain  rate, 

p  _3JL  i 

ijkJ  3S . .  k» 

x  '  7e - if  ,  jr  ;  JT  <49) 

<Eijkl  51^  H8Slj>  3Sij 

*  P 

Eliminating  A  in  the  expression  for  ekJl  in  Equation  4  8  then 
leads  to  the  rate  relation 


8ij  *  °ijk£  eki 


where  the  tensor  is  given  by 


Di jkl  “  Ei jki  "  BSijSkl 


and 


8 


0,  (elastic) 

3G/ (1+H/2G) k2  (plastic) 


(50) 


(51) 


(52) 


It  remains  to  determine  the  expansion  of  the  yield  surface 
as  plastic  flow  progresses  in  a  strain-hardening  material.  This 
information  is  obtained  by  requiring  that  any  point  in  stress 
space  remain  on  the  yield  surface  during  inelastic  deformation 
(that  is,  that  Equation  44  always  be  satisfied).  If  the  stress 
rate  s. .  are  known,  then 


e 

k 


1  3F  • 

2Tc  sij  * 


(53) 


In  practice,  the  elastic-plastic  constitutive  equation 
must  be  integrated  to  give  the  desired  relationship  between 
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stress  and  strain  increments  (see  Equation  22).  In  the  present 
formulation  this  integration  is  performed  using  a  trapezoidal 
rule  over  a  number  of  strain  subincrements  whose  size  is 
controlled  to  preserve  accuracy.  During  any  subincrement, 
then  Equations  50  and  53  become 

4,ij  ■  IEijk*  •  6*i'j  *ki>  4eki  154 

and 

Ak  *  5E  8ij  A8ij*  (55 

Depending  upon  the  size  of  the  computed  strain  increment. 
Equations  54  and  55  may  be  applied  tens  or  even  hundreds  of 
times  in  succession  to  update  the  element  stress  state  with  an 
acceptable  degree  of  accuracy. 


SECTION  3 

FINITE  ELEMENT  APPROXIMATIONS 

The  present  development  includes  four  classes  of  finite 
element  approximations:  thin  shell  or  sandwich  face  sheet- 
sandwich  core;  stiffening  members;  and  general,  three- 
dimensional  solids.  The  first  three  of  these  finite  element 
types  are  most  often  used  specifically  for  modeling  the 
individual  components  of  stiffened  sandwich  panels,  while  the 
solid  element  is  used  for  other  structure,  for  connections 
between  sandwich  and  other  components,  and  for  more  detailed 
analysis  of  sandwich  constructions  themselves. 

Each  of  the  finite  element  approximations  described  in 
this  Section  is  obtained  from  the  nonlinear  equilibrium  of 
Section  2  by  appropriate  specialization.  Linear  forms  of  the 
governing  equations  are  also  obtained  from  the  general 
formulation  through  the  assumptions  of  infinitesimal  displacements 
and  elastic  material  response. 

3.1  THIN  SHELL/FACE  SHEET  ELEMENTS 

The  formulation  of  Reference  11  is  used  herein  to  obtain 
a  thin  plate  and  shell  finite  element  suitable  for  representing 
sandwich  face  sheet  materials.  This  element  is  based  upon  a 
penalty  function  formulation,  which  leads  to  an  effective  thin 
shell  approximation  derived  directly  from  the  field  equations 
of  a  three-dimensional  continuum.  The  shell  element  so  derived 
possesses  several  distinct  advantages  over  similar  elements 
based  upon  a  specific  theory  of  thin  shells.  In  particular, 

-  geometric  parameters  of  an  element  are  completely 
defined  by  the  nodal  coordinate  information, 

-  complete  compatibility  of  displacements  with  standard 
isoparametric  solid  elements  is  possible, 

-  no  restrictions  need  be  placed  upon  the  extent  of 
displacements  or  rotations  of  the  shell,  and 
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-  effects  such  as  variable  thickness,  complex  shell 
intersections,  skewed  lateral  boundaries  and 
connections  to  adjacent  structure  are  accounted  for 
simply  and  effectively. 

The  undeformed  geometry  of  a  single  element  of  a  shell 
(or  face  sheet)  is  shown  in  Figure  2.  Geometric  parameters  are 
analyzed  with  respect  to  a  local  coordinate  system  (x,y)  which 
is  close  to  the  element  midsurface  z  =  z(x,y)  at  all  points. 

An  additional  coordinate  r,  is  used  to  describe  the  distance 
away  from  the  midsurface  at  any  point. 

It  is  assumed  that  on  the  element  level,  the  shell  is 
shallow,  so  that 


->  (56) 

z.y  1 

everywhere.  The  unit  vector  normal  to  the  shell  midsurface  at 
a  point  can  then  be  written  as 


n  = 


+  fc 


(57) 


in  the  local  system  of  coordinates.  Due  to  the  assumption  of 
shallowness,  the  (x,y,c)  coordinates  constitute,  at  least 
approximately,  a  local  Cartesian  systam;  taking  advantage  of 
this  fact,  the  position  vectors  of  an  arbitrary  point  within 
the  shell  before  and  after  deformation  may  be  expressed  in  the 
form 


(58) 


and 
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$  =  (x  -  CZ  „  +  u)i  +  (y  -  cz  v  +  v)j 
/X  » y 


+  (z  +  t,  +  w)k. 


(59) 


Equations  58  and  59  lead  directly  to  the  appropriate  definitions 
of  strain  for  the  shell  element,  since12 


ei;.dxidx;.  =  !l(d$*d$)  -  (3r*dr)]. 


(60) 


Thus,  the  Green *s  strains  referred  to  the  system  (x,y,C)  become: 


e  =  u  (1  -  ?z  )  -  Zz  v  +  z  w 
xx  9 xx  9 xy  f x  9 x  9 x 


+  i(u2  +  v2  +  w2  ) 

'  O  '  v  v  v' 


2'  ,x 


,x' 


e  =  v  11  -  tz  )  -  Cz  u  +zw 

yy  »yl  ,yy  »xy  *y  *y 


*  Ky  +  \y  +  w!y> 


ecc  w,c  “  Z,XU,C  "  z,yv,e 


+  l(u!c  +  +  w!t> 


2V  ’  \c(1  -  ”.yy’  +  \y  *  z.*u,y  +  z,y<w,E  '  v,y’ 


-  rz  u  +uu  +v  v_+w  w  _ 
^  ,xy  ,c  ,y  ,y  ,y  ,C 
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2e  *  u  (1  -  £z  )  +  v  (1  -  T,z  )  +  z  w  +  z  w 
xy  ,y  s  ,xx  ,x'  s  ,yy  ,x  ,y  ,y  ,x 


-  £Z  .  (u  +V„)+UU  +  V  V  +WW 

,xy  ,x  ,y  ,x  ,y  ,x  , y  ,x,y 


It  should  be  noted  that  the  strain-displacement  relations  above 
do  not  involve  the  orientation  of  the  displaced  midsurface 
normal,  and,  therefore,  apply  to  arbitrarily  large  displacements 
and  rotations  of  the  shell.  The  definitions  of  the  linear  and 
nonlinear  portions  of  the  incremental  strains  (A£. .  and  An^  •  in 
Equations  18  and  19)  follow  directly  from  Equations  62 

In  order  to  ensure  proper  behavior  of  the  finite  element 
approximation,  constraints  must  be  applied  to  enforce  the 
conditions  of  thin  shell  response.  To  prescribe  the  condition 
that  line  elements  initially  normal  to  the  shell  midsurface 
remain  straight,  normal  to  the  surface  and  unstrained,  it  is 
necessary  to  suppress  the  virtual  work  due  to  transverse  shear 
and  normal  strains,  and  to  enforce  constraint  relations  of  the 
form 


4ext  ‘  4eyC  =  4e„  ’  °-  (62> 

A  particularly  effective  method  of  enforcing  the  constraints 
of  Equation  62  is  developed  in  Reference  11,  by  defining  a  penalty 
functional  which,  in  the  limit,  forces  Equation  62  to  be  satisfied 
at  selected  discrete  points  within  an  element.  Define  the 
functional 
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4"c  ■ 


+  K4eyc<Jm> 


+  5VeCC(?mn 


in  which  M  is  the  number  of  constraint  locations,  W  ,W ,W  are 

positive  weights,  and  rm  is  the  location  at  which  the  m1-" 

constant  is  to  be  enforced.  Observe  that,  by  setting  6  (AIT  )  =  0, 

Equation  62  is  satisfied  identically  at  each  of  the  points 

r  ;  m=l,2 , . . . ,M. 
m 

Introducing  the  above  constraints  into  the  incremental 
principle  of  virtual  work,  Equation  21,  yields 


[D 

a 


BXu(4eX„  ‘  “x/T)44eaB  +  lSc.BS4r,aSldV 


-  At.6Au.dA  +  6  (All)  *  0. 

XX  C 


In  Equation  64,  Greek  subscripts  imply  a  range  of  summation  of  2, 
while  Latin  indices  have  a  range  of  3.  The  additional  thin 
shell  constraint  of  vanishing  normal  stress  through  the  thickness 
(o  *  0)  is  enforced  exactly  in  Equation  64.  Thus,  the  stress- 
strain  relation  for  the  remaining  nonzero  stresses  and  strains 
is  one  of  plane  stress. 

To  cast  Equation  64  in  matrix  form,  the  following 
quantities  are  defined.  Let 

uT  =  [_Au  Av  AwJ  (65) 

FT  *  Au  Au  Au  _  Av  „  Av  „  Av  r  Aw  „  Aw  Aw  _  I  (66) 
l»  rx  iY  »*>  i*  iY  tx  i y  » s J 

?T  -  L4exx  4eyy  24exy J  <67> 
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*  a 


a 


T 


W  * 


Liex{  4»»c 

L?xx  °yy  0xy-l 

0  0  ~ 

xx  xy  0 

axy  °yy  0 

L°  0  l 

“  x  0  0 

0  T  0 

0  0  T 

|_ot  a  aj 

"Wx  0  °“ 

o  w  o 

I  y 

jj)  0 


The  constrained 


tT  ■  I  At  At  At^  j 
~  •—  x  y  t,  — 1 

principle  of  virtual  work  then  becomes 
1 6eTD (e  -  aAT)  +  6FTx*F]dV 


M 


+  m=l  *1  (rm)  ~  yJr m)  " 


6uTxdA, 


3V 


(68) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 


(75) 


which  is  the  basis  of  the  thin  shell  finite  element.  It  remains 
only  to  specify  an  appropriate  approximation  for  u  in  terms  of 
a  finite  number  of  unknown  parameters. 
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In  selecting  the  finite  element  approximation  for  u  in 
Equation  75,  it  is  noted  that  only  first-order  derivatives 
of  the  displacement  appear  in  the  energy  form  of  the  equations. 
Thus,  any  approximation  which  preserves  continuity  of 
displacements  (but  which  may  violate  the  continuity  of 
slopes  between  elements)  is  admissible.  For  the  present 
development  the  simplest  possible  forms,  piecewise  liner 
polynomials,  are  selected.  Defining  the  natural  coordinates 
£i;i=»l,2,3  which  vary  between  -1  and  1  within  an  element,  the 
appropriate  one-dimensional  functions  are 


| (i  -  o 

Lo(0  =  i(i  +  n . 


Letting 

Nl(5)  * 

n2(§)  = 

n3(§)  = 
VP  = 
VP  = 

V§>  58 

N?  (p  = 


V?1> 

Ll(52) 

W 

w 

W 

W 

L2(^) 

l2(S2) 

L1 (^3 ) 

Ll(el> 

W 

Li(C3) 

Ll(?l) 

W 

W 

W 

W 

W 

W 

W 

W 

W 

W 

l2(?3) 

(76) 


(77) 


the  form  of  the  approximate  displacement  field  within  a  single 
element  is 
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U(C) 


(78) 


*  *^(5)  Ui  ■  NTU. 

Here  Ih;  i,2,...,8  are  values  of  the  displacement  u  at  the 
eight  nodal  points  of  an  element,  as  shown  in  Figure  3.  The 
displacements  u,  w  are  represented  in  exactly  the  same  manner. 

Substitution  of  the  above  approximation  into  Eguation  75 
leads  to  the  definition  of  the  tangent  stiffness  matrix  K^,  the 
geometric  stiffness  K  ,  a  discrete  "constraint  stiffness"  K  , 

-U  ''C 

and  the  external  force  vector  T.  In  matrix  form,  the 
equilibrium  relation  for  a  single  element  becomes 

(Kt  +  Kg  +  Kc)  U  =  T  (79) 

which  can  be  solved  iteratively  as  shown  in  Section  4. 

3.2  SANDWICH  CORE  ELEMENTS 

The  sandwich  core  finite  elements  employed  in  the  present 
work  are  obtained  directly  from  the  three-dimensional  nonlinear 
formulation  of  Section  2,  without  any  specialization.  Orthotropic 
stress-strain  relations,  as  described  in  Paragraph  2.2,  are 
permitted  for  analyzing  orthotropic  cores  (e.g. ,  metal  honeycomb 
core) ,  and  to  represent  the  very  low  extensional  moduli  typical 
of  many  types  of  core.  Permanent  deformations,  such  as  core 
crushing,  can  be  treated  in  an  approximate  fashion  using  the 
elastic-plastic  material  law  described  in  Paragraph  2.3. 

The  core  element  is  shown  in  Figure  4.  In  the  planform 
directions,  a  linear  interpolation  of  the  displacements  is  used 
for  compatibility  with  the  face  sheet  element  (Paragraph  3.1). 
Linear  displacement  shape  functions  are  also  employed  through  the 
core  thickness;  due  to  the  typical  high  flexibility  of  the  core 
layers  in  transverse  shear,  linear  displacements  are  nearly 
always  an  adequate  approximation.  Due  to  the  very  simple  nature 
of  the  displacement  state  approximation,  the  resulting  element 
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may  become  excessively  "stiff"  when  integrated  exactly**.  For 
this  reason,  a  single  integration  point  is  used  in  each  core 
element  to  provide  improved  accuracy  with  this  simple 
approximation.  Although  the  core  element  alone  would  prove  to 
be  unstable  using  the  single  point  integration  rule,  the 
addition  of  face  sheet  and  closeout  elements  provides  a 
stabilizing  effect  which  eliminates  singularity  in  the  final 
system  of  equations.  When  stabilized  in  this  manner,  the  simple 
core  element  used  here  provides  a  highly  effective  and  accurate 
representation  of  the  sandwich  response. 

3.3  STIFFENING  MEMBERS 

Beam  stiffener  elements  obtained  here  from  the  general 
formulation  of  Section  2  are  of  two  types:  face  sheet 
stiffeners  and  full-depth  stiffeners.  The  face  sheet  stiffening 
member  is  used  to  represent  discrete  beam  stiffeners  attached 
to  the  inner  or  outer  surface  of  a  sandwich  panel.  The  full- 
depth  stiffener  is  used  to  model  spar  panels  connected  to  both 
the  upper  and  lower  face  sheets  of  a  sandwich,  or  to  represent 
panel  edge  closeouts.  Both  of  these  configurations  are  shown 
in  Figure  5. 

Stiffener  elements  of  both  types  have  similar  geometry, 

shown  in  Figure  6.  The  coordinates  X^,  are  local  coordinates, 

embedded  in  the  plane  of  an  element.  In  this  local  coordinate 

system,  the  transverse  stresses  o  ,  a  and  a  are  assumed 

xz  yz  zz 

to  vanish,  and  thus  the  stress-strain  relation  for  an  element 
has  the  form  appropriate  for  a  plane  state  of  stress.  The 
resulting  element  models  extension,  primary  bending,  and  twisting 
effects.  Again,  compatibility  is  achieved  with  the  other 
sandwich  elements  through  the  use  of  linear  displacement  shape 
functions.  To  avoid  excessive  stiffness  of  the  element  in 
bending,  a  selective  integration**  scheme  is  used  for  the  in-plane 
shear  strain  energy. 


M 


FULL-DEPTH 
STIFFENER 


FACE  SHEET 
STIFFENER 


Figure  5.  Face  Sheet  Stiffener  and  Full-Depth  Stiffener 
Element  Configurations. 


3.4  THREE-DIMENSIONAL  SOLID  ELEMENTS 


A  three-dimensional  continuum  finite  element  is  included 
in  the  present  formulation,  for  use  in  the  modeling  of: 

-  general  structural  configurations, 

-  transitions  between  sandwich  panels  and  surrounding 
components,  and 

-  detailed  aspects  of  the  deformation  of  sandwich  panels. 

The  element,  which  has  a  variable  number  of  nodes,  can  be  made 
compatible  with  the  sandwich  face  sheet  and/or  core  elements 
(which  use  linear  shape  functions)  on  selected  element  surfaces 
while  higher-order  displacement  functions  are  used  elsewhere 
in  the  element.  Thus,  transitions  from  sandwich  panels  to 
attachments  or  other  supporting  structure  can  be  represented 
quite  simply.  Since  the  solid  element  is  fully  three-dimensional, 
it  is  also  useful  for  detailed  analysis  of  effects  such  as  core 
crushing  or  plastic  deformation  of  the  sandwich  face  sheets. 

The  general  solid  element  is  pictured  in  Figure  7.  Node 
points  1  through  8  (the  vertices)  are  required  for  each  element; 
however,  each  of  the  remaining  12  midside  nodes  may  be  included 
or  deleted  as  required.  In  most  instances,  the  use  of  the 
midside  nodes  is  highly  desirable  in  terms  of  accuracy.  The 
solid  element  may  be  numerically  integrated  using  either  an 
eight-point  (2x2x2)  or  27-point  (3x3x3)  Gaussian  quadrature 
formula. 
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SECTION  4 

NUMERICAL  SOLUTION  TECHNIQUES 

In  this  section,  the  methodology  used  in  solving  the 
governing  equations  for  a  complete  assemblage  of  finite  elements 
is  described.  Static  response  may  be  computed  with  or  without 
nonlinear  effects  included,  and  considering  both  thermal  and 
mechanical  loadings.  The  option  to  compute  a  specified  number 
of  natural  vibration  frequencies  and  their  associated  mode 
shapes  is  also  provided. 

4.1  STATIC  RESPONSE  SOLUTION 

In  static  problems,  the  finite  element  equations  take 
the  form 

(Kt  +  Kg)AU  =  AT  (80) 

in  which 

Kt  =  tangent  stiffness  matrix, 

Kq  =  geometric  (initial  stress)  stiffness, 

AU  =  incremental  nodal  displacements,  and 

AT  =  incremental  nodal  forces,  including 
both  mechanical  and  thermal  effects. 

The  tangent  stiffness,  for  nonlinear  problems,  is  a  function  of 
the  total  displacements,  KT  +  KT(U),  and  the  geometric  stiffness 
depends  upon  the  current  stress  level,  =  KG  ( a ) . 

When  the  problem  is  assumed  to  be  linear  (small  displacements, 
elastic  material  behavior),  K„  =  K  ,  the  linear  stiffness,  and 

—  1  ~  £j 

KG  =  0.  In  this  case,  Equation  80  can  be  written  using  the  total 
displacements  U  and  forces  T: 

Ke  U  =  T.  (81) 
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A  direct  solution  is  obtained  by  factoring  K£  as 

Ke  =  LDLT  (82) 

where  L  is  a  unit  lower  triangular  matrix  and  D  is  diagonal. 

Next  letting 

Z  =  DLTU,  (83) 

the  relation 

LZ  =  T  (84) 

can  be  solved  directly  for  Z.  Similarly,  Equation  83  gives 

LTU  =  D_1Z  (85) 

which  can  finally  be  solved  for  U.  Knowing  the  displacement 
solution,  strain  and  stress  information  can  be  computed  at 
selected  locations  in  the  model  to  complete  the  linear  analysis. 

When  nonlinear  effects  are  considered  in  the  solution 
of  Equation  80,  an  iterative  technique  must  be  adopted.  In 
addition,  the  nonlinear  solution  is  usually  performed  in  an 
incremental  fashion,  by  applying  the  total  loading  in  a  number 
of  steps.  Such  an  approach  is  adopted  to  provide  the  static 
solution  for  an  entire  range  of  loading,  to  enhance  convergence, 
and  to  facilitate  the  detection  of  buckling  or  similar 
phenomena.  In  plastic  analysis,  this  incremental  method  is 
necessary  due  to  the  history-dependent  nature  of  the  response. 

The  nonlinear  equations  are  solved  using  various  forms  of 
the  Newton-Raphson  method,  as  shown  in  Figure  8.  First,  the 
nonlinear  system  is  written  as 

(Kt  +  Kg)AU  =  T  -  I  (86) 
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MODIFIED  NEWTON- RAPHSON 
(CONSTANT  STIFFNESS) 
ITERATION 


DEFLECTION 


FULL  NEWTON -RAPHSON 
ITERATION 


X, 


where  T  is  the  total  applied  loading  and  I  are  the  internal 
forces  based  upon  the  current  (computed)  state  of  stress.  Once 
the  computed  stress  state  is  in  equilibrium  with  the  external 
loading,  the  right-hand  side  of  Equation  86  becomes  zero.  Thus, 
Equation  86  can  be  used  repeatedly  to  compute  corrections  to  an 
estimated  displacement  state  for  given  T,  until  the  equilibrium 
condition  (T  »  I)  is  satisfied  to  within  a  specified  tolerance. 
During  a  full  Newton -Raphson  solution,  the  coefficient  matrix 
is  updated  at  each  iteration  cycle.  In  the  modified  Newton 
solution,  the  stiffness  matrix  is  reformed  only  at  the  start 
of  an  increment,  as  indicated  in  Figure  8.  The  conventional 
Newton's  method  tends  to  converge  rapidly,  with  each  iteration 
consuming  a  relatively  large  amount  of  processor  time;  with  the 
modified  iteration,  each  cycle  of  iteration  is  computationally 
simpler,  but  convergence  is  generally  rather  slow. 


4.2  NATURAL  FREQUENCY  SOLUTION 

When  free  vibrations  are  considered,  the  problem  is  assumed 
to  be  linear,  and  the  forcing  function  T  is  replaced  by  the 
inertial  forces, 

K£U  =  -MU  (87) 

where  M  is  the  structure  mass  matrix.  For  harmonic  motions, 
then,  the  nodal  displacements  are 


U  =  X  sinu)T 

and  the  equation  of  motion  becomes 

KeX  =  U)2MX. 


(88) 


(89) 


Notice  that  the  vector  X  describes  only  the  relative  displacements 
of  the  node  points,  and  represents  a  mode  shape  corresponding  to 
the  circular  frequency  of  vibration,  u>.  If  K  and  M  are  matrices 
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of  order  N,  then  there  exist  N-solutions  of  Equation  89:  uk  , 

Xi;  i=l,2...,N. 

Generally  only  the  lowest  few  frequencies  and  mode 

shapes  of  the  system  will  be  of  interest,  and,  therefore,  a 

vector  iteration  procedure  is  most  aopropriate.  In  the  present 

13 

development,  the  method  of  simultaneous  vector  iteration  is 
used. 

If  KE  can  be  expressed  in  the  Choleski-factored  form 

KE  =  LLT  (90) 


with  I,  a  lower  triangular  matrix,  then  the  transformation 


Y 


T 

LX 


(91) 


allows  Equation  89  to  be  written  in  the  standard  form 


(92) 


Here  A  is  defined  by 


A  =  L_1ML_T.  (93) 

Given  a  set  of  trial  solution  vectors  Y.;  i=l,2,...,m, 
a  partial  modal  matrix  is  formed, 

*  -  [Y.  Y,  •  •  •  Y  ]  .  (94) 

~  £  ~m 

Noting  that  any  set  of  eigenvectors  of  A  must  be  orthogonal 
with  respect  to  A,  the  next  step  is  to  inspect  the  "interaction 
matrix" 

B  =  $TA4>  (95) 
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which  is  diagonal  if  the  columns  of  $  are  true  eigenvectors. 

When  the  vectors  are  properly  normalized,  the  diagonal 
entries  of  B  are  estimates  of  the  eigenvalues 

uk  -  /bii  (96) 

Before  the  solution  has  converged,  the  interaction  matrix 
B  is  not  diagonal.  Instead,  an  approximate  solution  can  be 
performed  for  B  to  obtain  an  improved  estimate  of  the  eigenvalues. 
The  trial  vectors  (i.e.,  columns  of  ♦)  are  then  modified 
accordingly,  orthogonalized,  and  used  in  the  next  iteration 
cycle. 

A  solution  for  the  lowest  several  frequencies  and  mode 
shapes  using  simultaneous  iteration  tends  to  converge  in  very 
few  iterations,  provided  the  number  of  trial  iteration  vectors 
is  sufficient.  If  p  is  the  number  of  natural  frequencies  to 
be  solved,  an  estimate  for  the  proper  number  of  trial  vectors 
can  be  obtained  from 

m  =  minfZp,  p  +  5] .  (97) 

The  procedure  in  general  tends  to  be  quite  fast,  since  only  one 
factorization  of  KE  is  required,  and  can  be  performed 
economically  without  actually  forming  the  matrix  A  (Equation  92). 


SECTION  5 
SAMPLE  ANALYSES 


Several  sample  problems  are  presented  in  this  Section  to 
demonstrate  the  analytical  capabilities  of  the  present  finite 
element  formulation.  Classes  of  problems  considered  include 
large  deflection  and  elastic-plastic  analysis,  buckling,  natural 
frequency  solutions  and  thermal  stress  analysis. 

5.1  NONLINEAR  ANALYSIS  OF  A  SANDWICH  PANEL  UNDER  PRESSURE 

A  square  sandwich  panel,  50  inches  on  each  side,  is 
subjected  to  a  uniform  lateral  pressure.  The  three-layer  plate, 
shown  in  Figure  9,  has  identical  aluminum  face  sheets 
(E  =  10.5  x  10^  psi;  v  =  0.3)  0.015  inches  thick,  bonded  to 
an  aluminum  honeycomb  core  of  one  inch  thick.  The  core  is 
assumed  to  be  isotropic,  with  shear  modulus  G  =  50,000  psi. 

All  boundaries  of  the  sandwich  are  fully  clamped. 

Due  to  symmetry  of  the  geometry  and  loading,  one  quarter 
of  the  panel  is  considered  in  the  finite  element  solution. 

The  finite  element  discretization  consists  of  a  total  of  75 
finite  elements,  25  in  each  layer.  The  two  face  sheets  are 
modeled  using  eight-node,  thin  shell  elements  (Paragraph  3.1). 
Three-dimensional  sandwich  core  elements  (Paragraph  3.2)  are 
used  for  the  central  layer.  Note  that  these  element  types  are 
fully  compatible  so  that  no  special  constraints  are  necessary 
for  joining  the  individual  layers.  The  nonlinear  solution 
has  been  obtained  in  load  increments  of  one  psi  to  a  total 
pressure  of  20  psi  followed  by  two  psi  increments  to  30  psi. 

The  nonlinear  central  displacement  of  the  sandwich  is 

plotted  versus  load  in  Figure  10.  Nonlinear  finite  element 

2 

results  obtained  by  Monforton  ,  using  16  specially  formulated 

bicubic  sandwich  elements,  are  shown  for  comparison.  Agreement 

between  the  two  finite  element  solutions  is  quite  good.  Figure 

14 

10  also  shows  the  perturbation  solution  of  Kan  and  Huang  , 
given  by 
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Sandwich  Panel  Under  Uniform  Lateral  Pressure. 


(98) 


q  =  10.5299W  +  4.8550w3 

c  c 

in  which  q  is  the  applied  pressure  and  wc  the  transverse  center 
displacement.  The  analytical  solution  of  Reference  14  is  valid 
for  deflections  which  are  smaller  than  the  core  thickness,  and 
reasonable  agreement  with  the  two  numerical  solutions  is  observed 
in  this  region.  For  larger  deflections,  the  perturbation 
analysis  requires  more  terms  for  acceptable  accuracy;  the  two- 
term  solution  gives  results  which  overestimate  the  influence  of 
membrane  stiffening  upon  the  panel  deflection. 

5.2  NATURAL  FREQUENCIES  OF  SANDWICH  PANELS  WITH  VARIOUS 
BOUNDARY  CONDITIONS 

The  three-layer  sandwich  plate  shown  in  Figure  11  has 
been  analyzed  for  its  free  vibration  response  using  a  number  of 
different  edge  conditions.  The  finite  element  model  of  the  panel 
consists  of  36  elements  in  each  layer,  and  contains  a  total  of 
196  nodes  and  480  unconstrained  degrees  of  freedom.  Material 
properties  and  dimensions  of  the  specimen  are  as  follows: 


Dimensions  : 

a 

S3 

62.25 

in. 

b 

= 

43.50 

in. 

*f 

= 

0.072 

in. 

fcc 

= 

1.856 

in. 

Face 

Sheets : 

E 

= 

1.0  x 

107 

V 

= 

0.  33 

Core 

Layer  : 

Gxz 

= 

30000. 

psi 

Gyz 

= 

30000. 

psi 

Ez 

= 

60000. 

psi 

Three  types  of  edge  conditions  have  been  considered  in  the 
present  analysis.  The  first  is  a  simply-supported  condition  in 
which  in-plane  motions  of  the  edqe  are  permitted;  this  type  of 
constraint  allows  a  pure-bending  type  of  response  of  the  entire 
panel  and  is  readily  verified  with  known  analytical  solutions. 
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For  the  second  type  of  constraint,  in-plane  displacements  are 
prevented  on  the  lower  face  of  the  panel ,  introducing  a  net 
stretching  effect  when  transverse  bending  occurs.  The  third 
boundary  condition  used  is  a  fully  clamped  edge. 

Fundamental  frequencies  of  these  panels  are  summarized 
in  Table  1.  For  the  first  set  of  edge  conditions,  analytical 
frequencies  computed  from  the  formulation  of  Reference  15  are  also 
listed.  Agreement  between  the  two  solutions  is  good.  The 
stiffening  effect  produced  by  the  two  alternative  boundary 
conditions  is  clearly  shown  in  the  frequency  results. 

5.3  BUCKLING  OF  A  SIMPLY-SUPPORTED  SANDWICH  PLATE 

The  compressive  buckling  of  a  plate  of  sandwich 
construction  is  considered.  A  square,  three-layer  panel  (Figure 
12)  is  subjected  to  a  uniform  compressive  load  of  Nx  pounds  per 
inch.  The  panel  is  23.5  inches  on  each  side,  and  supported  at 
each  face  of  the  sandwich  on  all  four  edges  (vertical 
displacements  only  are  prevented) . 

The  outer  face  sheets  of  the  panel,  which  are  represented 
by  thin  shell  finite  elements,  are  each  0.021  inches  in  thickness, 
with  isotropic  material  properties  E  =  9.5  x  106  psi,  v  =  0.3. 

The  core  layer,  0.181  inches  thick,  has  a  transverse  shear 
rigidity  G  =  19000  psi.  Each  layer  of  the  model  contains  16 
elements  of  equal  planform  dimensions.  Only  one  quadrant  of 
the  panel  is  considered  in  the  numerical  solution,  due  to 
symmetry  of  the  geometry  and  loading.  On  the  lateral  boundaries, 
the  tangential  transverse  shear  strains  within  the  core  are 
suppressed  by  making  the  upper  and  lower  face  sheet  displacements 
equal  in  the  direction  parallel  to  each  edge. 

A  solution  for  the  buckling  load  NCR  has  been  obtained 
by  applying  the  in-plane  forces  incrementally  until  a  sudden 
increase  in  transverse  displacement  is  observed.  Out  of  plane 
deflections  are  triggered  by  a  small  (one  pound)  transverse 
load  applied  at  the  center  of  the  plate.  Buckling  is  found  to 
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TABLE  1 


NATURAL  FREQUENCY  RESULTS  FOR  FLAT  SANDWICH  PANEL 
WITH  VARIOUS  BOUNDARY  CONDITIONS 


Boundary  Condition 

Frequency,  Hz 
(Computed) 

. 

Frequency,  Hz 
(Analytical) 

Simple  Support 

Inplane  Motion  Permitted 

174 

170 

Simple  Support 

No  Inplane  Motion 

211 

- 

Clamped  Edge 

288 

i 
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occur  for  an  applied  load  of  =  305  pounds  per  inch;  this 
computed  value  compares  well  with  previous  analytical  and 
experimental  results,  as  shown  in  Table  2.  It  is  noted  that  all 
of  the  analytical  results  give  estimates  of  the  critical  load 
which  are  about  nine  percent  too  high;  it  is  likely  that  the 
assumption  of  zero  transverse  shear  strains  at  the  panel 
boundaries  is  largely  responsible  for  this  error.  The 
transverse  shear  constraint  has  been  used  in  the  present  finite 
element  solution  to  permit  comparison  with  previous  analytical 
solutions. 

5.4  PLASTIC  ANALYSIS  OF  A  TOROIDAL  SANDWICH  SHELL 

A  toroidal  sandwich  shell  panel  under  line  loading  has  been 
analyzed  to  determine  its  elastic-plastic  response.  The  following 
geometric  parameters  are  used: 

major  radius  R  =  5.0  inches 

included  angle  0O  =  75.0  degrees 

minor  radius  r  =  1.5  inches 

included  angle  9r  =  90.0  degrees 

core  thickness  t  =0.90  inches 

c 

face  thickness  t^  =  0.05  inches 

The  geometry  of  the  shell  is  shown  (in  the  form  of  a  finite 
element  model)  in  Figure  13.  Figure  14  shows  the  outer  face 
sheet  of  the  sandwich  with  node  numbers  on  the  exterior  surface 
labeled.  A  line  load  is  applied  in  the  vertical  direction  along 
the  upper  edge  of  the  shell  (nodes  36,  56,  92,  ...,  316  in 
Figure  14).  The  two  short  edges  of  the  panel  are  fully  clamped. 

The  panel  is  constructed  of  aluminum,  with  the  properties 
Ef  =  1.0  x  107  lb/in2 
Vf  =  0.30 

Gq  =  30000.  lb/in2 
2 

A  yield  stress  of  10000  lb/in  is  assumed  in  all  layers,  and 
the  material  is  considered  to  be  elastic,  perfectly  plastic. 
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TABLE  2 

COMPARISON  OF  BUCKLING  LOADS  FOR  SIMPLY-SUPPORTED  SANDWICH  PANEL 


Reference 

Method 

Ncr  (pounds/inch) 

Hoff16 

Series  Solution 

303.0 

Plantema16 

Series  Solution 

308.0 

17 

Brockman 

Series  Solution 

309.0 

Boiler16 

Experimental 

266.  -  300. 

2 

Mon f orton 

Finite  Element 

307.5 

Present  Analysis 

Finite  Element 

305.0 
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pper  Face  Sheet  of  Toroidal  Shell  Panel  with  Node  Numbers 


From  an  elastic  analysis  of  the  shell,  the  yield  stress 
is  found  to  occur  at  a  total  load  of  approximately  600  pounds. 

The  elastic-plastic  solution  is  then  performed  using  a  single 
increment  of  600  pounds,  followed  by  equal  increments  of  180 
pounds  to  a  maximum  loading  of  1500  pounds.  Full  Newton- 
Raphson  iterations  are  used  at  each  step  of  the  solution  to 
maintain  equilibrium  of  the  nonlinear  system. 

The  transverse  displacement  at  the  center  of  the  loaded 
edge  (node  168  in  Figure  14)  is  tabulated  as  a  function  of  the 
total  load  in  Table  3.  Deflected  shapes  of  the  inner  and 
outer  face  sheet  layers  are  also  given  in  Figures  15  and  16, 
respectively,  for  the  maximum  value  of  loading.  Although  the 
nonlinearity  in  central  displacement  at  this  loading  level 
are  rather  mild  (see  Table  3) ,  material  yielding  is  extensive 
near  the  clamped  edges  and  the  loaded  free  boundary.  Figure  17 
shows  the  zone  of  plastic  behavior  in  the  outer  face  sheet  at 
the  maximum  loading  level. 

5.5  THERMAL  STRESS  ANALYSIS  OF  A  VARIABLE-THICKNESS  SANDWICH 

STRIP 

A  narrow  sandwich  strip  with  a  small  geometric  asymmetry 
is  subjected  to  uniform  heating,  to  study  the  effect  of  the 
asymmetry  in  producing  transverse  deflections  at  high  temperature. 
One  quarter  of  the  symmetric  panel  is  shown  in  Figure  18.  The 
uniform  middle  segment  of  the  strip  is  16  inches  in  length,  while 
the  tapered  sections  and  ends  are  each  two  inches  long.  The 
lower  face  sheet  is  0.075  inches  in  thickness.  At  the  small  ends, 
the  thicknesses  of  the  remaining  layers  are  t  =  0.25  inches, 
tj  =  0.050  inches,  increasing  to  tc  =  0.50  inches,  t^  =  0.075 
inches  in  the  central  section.  Full-depth  stiffeners  0.015  in 
width  are  used  to  close  out  the  lateral  boundaries  of  the  strip. 

The  entire  strip  is  constructed  from  aluminum,  the 
properties  being 


TABLE  3 


LOAD-DEFLECTION  HISTORY  AT  CENTRAL  POINT  ON 
LOADED  EDGE  OF  TOROIDAL  SANDWICH  PANEL 


16.  Deformed  Shape  of  Outer  Face  Sheet  at 
Maximum  Loading. 


Figure  17.  Plastic  Boundary  in  Outer  Face  Sheet  of  Toroidal 
Shell  at  Maximum  Load 


Face  Sheets: 


E 


10.  x  106  lb/in2 


Stiffeners : 


Core : 


v  =  0.30 

a  =  12.  x  10  6  in/in°F 

E  =  10.  x  106  lb/in2 
v  =  0.30 

a  =  12.4  x  10  6  in/in°F 

G  =  20000.  lb/in2 
X  z 

G  =  30000.  lb/in2 
yz 

ax=«  =  1.  x  10-10  in/in°F 

a  =  12.  x  10“®  in/in°F 
z 


The  finite  element  model  of  the  strip  is  shown  in  more 
detail  in  Figures  19  and  20.  Figure  19  shows  the  face  sheet 
and  full-depth  stiffener  elements;  Figure  20  is  a  plot  of  the 
core  layer  elements  only,  in  a  similar  orientation.  The  model 
represents  only  one  fourth  of  the  strip  due  to  double  symmetry 
of  the  problem.  The  small  ends  of  the  strip  are  fully  clamped 
at  all  points. 

A  deformed  geometry  plot  of  the  panel  (Figure  21)  after 
uniform  heating  shows  quite  clearly  the  bending  response 
caused  by  the  asymmetry  at  the  panel  ends.  The  displaced 
shape  (solid  lines)  is  drawn  to  scale  for  a  temperature  rise 
of  2000°F.  While  the  deflections  remain  linear,  the  transverse 

central  displacement  due  to  uniform  temperature  rise  is 

-4 

approximately  5.0  x  10  in/°F. 


SECTION  6 

SUMMARY  AND  CONCLUSIONS 

A  finite  element  approach  for  the  linear  and  nonlinear 
analysis  of  structures  incorporating  sandwich  composite  materials 
has  been  described.  The  methodology  is  based  upon  a  philosophy 
of  modeling  each  layer  of  a  sandwich  material  in  some  detail, 
rather  than  using  approximate  properties  assembled  from  the 
characteristics  of  each  layer.  With  such  an  approach,  it  is  a 
rather  straightforward  task  to  simulate  response  problems 
involving  local  deformations,  crushing,  buckling  and  plasticity. 
Furthermore,  connections  between  sandwich  panels  and  other  more 
conventional  construction  materials  may  be  represented  in 
detail  using  standard  continuum  finite  elements. 

The  present  formulation  includes  lower-order  thin  shell 
and  orthotropic  solid  elements  for  representing  the  basic 
characteristics  of  sandwich  panels,  and  compatible  bending 
elements  for  modeling  face  sheet  or  full  depth  stiffeners. 
Standard  three-dimensional  solid  elements  describe  other 
connected  structure,  or  may  be  used  for  more  detailed  analysis 
of  sandwich  panels  when  required.  Each  of  these  elements  is 
applicable  to  problems  involving  arbitrarily  large  displacements 
or  rotations,  since  no  restrictive  assumptions  are  made  in  the 
element  theoretical  formulations.  Plasticity  is  considered 
using  the  von  Mises  yield  condition  and  the  Prandtl-Reuss 
equations  of  incremental  plastic  flow. 

The  finite  element  approach  described  has  been  shown  to 
be  effective  for  numerous  applications,  involving  both 
mechanical  and  thermal  loading.  Natural  frequency  calculations 
indicate  that  the  methodology  is  also  quite  accurate  for  this 
purpose. 
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APPENDIX 

COMPUTER  PROGRAM  INPUT  INSTRUCTIONS 

The  theoretical  and  numerical  formulation  described  in  this 
report  have  been  implemented  in  a  computer  program,  from  which 
representative  solutions  have  been  presented  in  Section  5.  In  the 
following,  input  data  instructions  for  use  of  the  program  are 
provided. 

The  computer  program  input  is  divided  into  26  input  Data 
Sets,  not  all  of  which  will  be  required  in  a  single  analysis. 
Conditions  under  which  each  Input  Set  should  be  included  or 
skipped  are  noted  at  the  start  of  the  Data  Set.  Default 
values  and/or  additional  explanatory  information  concerning 
each  item  of  input  are  included  in  each  Input  Set  description 
in  the  form  of  Notes. 

For  each  item  of  input,  a  corresponding  FORTRAN  variable 
name  is  listed  which  serves  to  define  the  data  type  (floating 
point,  integer).  Any  exceptions  to  this  scheme  (e.g., 
alphanumeric  data)  are  noted  explicitly  as  they  occur.  All 
floating  point  (real)  data  may  be  entered  with  or  without 
exponents;  exponents,  if  used,  must  be  right- justified  in  the 
data  field  provided.  Integer  value  must  be  input  without 
decimal  points  or  exponents,  and  must  be  right- justified  in 
the  data  field. 

No  plotting  utilities  are  provided  in  the  present  program 
for  display  and  verification  of  the  finite  element  model  data. 

For  this  reason,  the  nodal  coordinate  and  element  connectivity 
data  have  been  made  compatible  with  that  used  in  the  MAGNA 
finite  element  program.  Therefore,  full  advantage  can  be  taken 
of  the  capabilities  offered  in  the  MAGNA  program  for  interactively 
displaying  the  model  geometry  for  verification  and  documentation 
purposes. 
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INPUT  SET  1 


PROBLEM  IDENTIFICATION 
(Required  for  all  analyses) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-80 

TITLE 

Alphanumeric  Problem  Description 

- 
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INPUT  SET  2 


SOLUTION  OPTIONS 
(Required  for  all  analyses) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

B 

IATYPE 

Analysis  Type 
=1,  Linear  Static 
=2,  Nonlinear  Static 
=3,  Natural  Frequency 

- 

6-10 

ITHERM 

Thermal  Stress  Flag 

=0,  Thermal  Stresses  Neglected 

=1,  Thermal  Stresses  Included 

11-15 

IP0ST 

Postprocessor  File  Flag 
=0,  Do  Not  Save  Results 
=1,  Save  Analysis  Results 
on  File  PLOTFIL 

2 

H 

NINCR 

Number  of  Solution  Increments 

- 

■ 

IPRNTF 

Printing  Frequency 

- 

3 

1-10 

STEP 

Loading  Parameter  Step  Size 

(1) 

NOTES : 

(1)  In  a  nonlinear  analysis,  all  external  forces 
are  expressed  as  functions  of  a  "loading 
parameters,"  which  is  continuously  increasing 
during  the  solution.  STEP  defines  the  amount 
by  which  the  loading  parameter  (whose  initial 
value  is  zero)  is  increased  at  each  solution 
increment. 
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INPUT  SET  3 


EQUILIBRIUM  ITERATION  OPTIONS 

(Required  only  for  nonlinear  analysis,  IATYPE  =  2  on  input  Set  2) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

ITTYPE 

Iteration  Type 
=0,  No  Iteration 
=1,  Modified  Newton  Iteration 
=2,  Newton-Raphson  Iteration 

(1) 

6-10 

ITFREQ 

Iteration  Frequency 

(2) 

11-15 

ITMAX 

Maximum  Number  of  Iterations/ 
Increment 

(3) 

16-25 

EQTOL 

Iteration  Convergence  Tolerance 

(4) 

NOTES i 

(1)  ITTYPE  determines  the  type  of  equilibrium  iteration 
(if  any)  to  be  performed  at  specified  intervals 

in  the  solution  to  restore  the  nonlinear  conditions 
of  equilibrium  (i.e.  ,  internal  forces  »  external 
forces) .  Modified  Newton -Raph son  iteration  (ITTYPE  = 

1)  involves  no  reformulation  of  the  system  stiffness 
matrix;  internal  forces  are  computed  at  each  iteration, 
and  the  resulting  out-of-balance  forces  are  applied 

to  obtain  displacement  corrections  to  improve  the 
solution.  With  full  Newton-Raphson  iteration  (ITTYPE  = 

2) ,  the  stiffness  is  formed  and  solved  at  each  iteration. 
This  procedure  is,  therefore,  more  expensive  than  the 
modified  iteration,  but  is  normally  quicker  to  converge. 

(2)  Iteration  is  performed  every  ITFREQ  increments  of  the 
nonlinear  analysis.  The  default  value  is  ITFREQ  =  1. 

(3)  If  more  than  ITMAX  iterations  are  required  in  any 
single  solution  increment,  the  solution  will  be 
terminated  to  permit  a  change  in  increment  size  or 
other  parameters.  The  default  value  is  ITMAX  =  20. 


I  / 
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(4)  EQT0L  defines  the  convergence  tolerance  on  errors  in 
the  internal  forces  R,  as  a  fraction  of  the  applied 
load,  F,  during  equilibrium  iteration.  The 
iteration  is  considered  converged  when 


II  p  -  R I  i 

- = - = -  <  EQT0L  , 

II  F  II 

where  ||  V  jj  denotes  the  Euclidean  norm,  /  VTV. 


INPUT  SET  4 


£ 


EIGENVALUE  SOLUTION  PARAMETERS 

(Required  only  for  natural  frequency  analysis,  IATYPE  =  3  in 
Card  Set  2) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

NTRIAL 

Number  of  Iteration  Trial 

Vectors 

(1) 

6-10 

NREQD 

Number  of  Natural  Frequencies 
to  be  Determined 

(2) 

11-15 

MAXIT 

Maximum  Number  of  Iterations 

(3) 

16-25 

T0LVEC 

Vector  Tolerance  for 

Convergence  of  Frequency 
Solution 

(4) 

NOTES : 

(1)  In  general,  the  greater  the  number  of  trial  iteration 
vectors,  NTRIAL,  the  better  the  convergence  and 
accuracy  characteristics  of  the  solution.  However, 
the  use  of  an  excessive  number  of  iteration  vectors 
is  costly  and  inefficient  in  terms  of  central  memory 
requirements.  The  use  of  the  NTRIAL  =  min  (2*N,N+5), 
where  N  is  the  number  of  frequencies  to  be  solved, 
provides  a  good  balance  between  rate  of  convergence 
and  storage  requirements.  The  default  value  is 
NTRIAL  =  2. 

(2)  The  program  will  determine  the  first  NREQD  natural 
frequencies  and  normal  modes  of  the  linear  system 
KX  =  u>2  MX,  where  K  is  the  system  stiffness  matrix 

and  M  the  mass  matrix.  Since  the  solution  is  performed 
by  vector  iteration,  NREQD  is  limited  to  values  which 
are  relatively  small  for  large  finite  element  models. 
The  default  value  is  NREQD  =  1. 

(3)  MAXIT  controls  the  total  number  of  iteration  cycles 
performed  during  the  solution.  A  value  of  MAXIT  =  15 
to  20  is  sufficient  for  nearly  all  problems,  unless  the 
number  of  frequencies  to  be  computed  (NREQD)  is  quite 
large.  A  default  value  of  MAXIT  =  15  is  used. 
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(4)  TOLVEC  defines  the  convergence  tolerance  on 

successive  approximations  to  each  eigenvector. 

If  and  are  successive  iterates  to  a  single 

eigenvector,  that  eigenvector  is  considered 
converged  if 

II  Vi+1  -  V.  ||  <  TOLVEC  , 

where  |jv||  denotes  the  Euclidean  norm,  /  yTy  .  The 
default'is  TOLVEC  *  0.001. 
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INPUT  SET  5 


1 

NODAL  COORDINATES 
(Required  for  all  analyses) 

The  number  of  cards  entered  in  this  Section  is  determined 
by  the  number  of  nodal  points  to  be  defined  in  the  model.  Nodal 
input  is  terminated  by  a  single  blank  card  (i.e.,  NODE  =  0). 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

Literal  "COORDINATES" 

- 

16-20 

NODES 

Total  Number  of  Node  Points 

- 

21-25 

(blank) 

26-35 

TDFLT 

Default  Nodal  Temperature 

Value 

(1) 

36-40 

ITDATA 

Load  Parameter  Curve  for 
Variation  of  Nodal 
Temperatures 

(2) 

2-n 

1-5 

NODE 

Node  Point  Number 

(3) 

6 

ISYS 

Reference  Coordinate  System 
=  :  Cartesian  X,Y,Z, 

=A:  Cylindrical  R,9,Z 
=B:  Spherical  R,<J>,0 

(4) 

7-10 

NINCR 

Increment  for  Node  Point 
Generation 

(5) 

11-20 

X ( NODE ) 

Coordinate  X^ 

21-30 

Y (NODE) 

Coordinate  X2 

31-40 

Z  (NODE) 

Coordinate  X^ 

41-50 

T (NODE) 

Nodal  Temperature 

(6) 
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NOTES : 


(1)  TDFLT  is  a  default  temperature  value  which  will  be 
assigned  to  any  node  point  whose  input  temperature 
value  is  zero.  This  parameter  is  commonly  used  in 
problems  of  uniform  heating,  etc.  ,  in  which 
temperature  data  may  be  absent  from  the  original 
input  data. 

(2)  ITDATA  refers  to  a  loading  parameter  curve  (i.e., 
function  of  time)  which  is  input  in  input  set  25. 

Curve  ITDATA  describes  the  variation  of  temperature 
at  all  nodes  in  the  model,  and  is  used  in  nonlinear 
analysis  only. 

(3)  Acceptable  nodal  point  numbers  are  between  1  and  NODES. 
Not  every  node  need  be  connected  to  an  active  element 
in  the  model,  but  inactive  nodes  must  be  fully 
constrained  (through  boundary  condition  input) . 
Coordinate  data  is  read  until  a  blank  (i.e., 

NODE  =0)  is  encountered. 

(4)  Nodes  may  be  defined  in  circular  cylindrical 
coordinates  by  setting  ISYS  =  A  and  providing  as 
input  the  R,0,Z  coordinates  of  the  point,  where  6 
is  measured  in  degrees.  In  this  case,  the  node 
coordinates  are  converted  internally  to  Cartesian 
coordinates  defined  by: 

X  =  Rcos9 

Y  =  Rsin9 

Z  =  Z. 

When  ISYS  =  B,  the  program  interprets  coordinate 
data  as  spherical  coordinate  values  R,  <J> ,  and  6, 
where  both  0  and  0  are  measured  in  degrees. 

Spherical  coordinates  are  then  converted  to  Cartesian 
coordinates  by  the  formulas: 

X  =  Rsin0cos0 

Y  =  Rsin0sin9 

Z  =  RCOS0. 

(5)  Node  generation  increments  NINCR  are  entered  on  the 
second  card  of  a  pair,  causing  nodes  to  be  equally 
spaced  between  the  last  and  current  nodes,  with 
numbering  increment  NINCR.  As  an  example,  the  data 
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10 


0. 


0. 


0. 


20  2  10.  -10.  0. 

is  equivalent  to 

10  0.  0.  0. 

12  2.  -2.  0. 

14  4.  -4.  0. 

16  6.  -6.  0. 

18  8.  -8.  0. 

20  10.  -10.  0. 

Note  that  incremental  node  generation  is  performed 
in  the  Cartesian  system  only. 

(6)  Nodal  temperatures  are  understood  to  be  the 

differences  in  temperature  from  the  (unstressed) 
reference  of  the  structure  (usually  "room  temperature"). 
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INPUT  SET  5 

SIZING  PARAMETERS  FOR  SANDWICH  CORE  ELEMENTS 

(Skip  this  Set  if  no  Sandwich  Core  Elements  are  to  be  defined) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

IECOPE 

Element  Type  Code?  Enter  the 
number  "2" 

6-10 

NMAT 

Number  of  Material  Property 

Sets 

(1) 

11-15 

NELEM 

Number  of  Elements  of  +-his 
Element  Type 

- 

16-20 

NAXIS 

Number  of  Orthotropic  Axis 
Definitions 

(2) 

NOTES : 

(1)  NMAT  defines  the  total  number  of  property  sets  to 
be  defined  in  Input  Sets  6  and  7. 

(2)  NAXIS  determines  the  number  of  axis  sets  to  be 
defined  in  Input  Set  8. 


INPUT  SET  6 


ISOTROPIC  SANDWICH  CORE  PROPERTIES 

(Skip  this  Set  if  no  Sandwich  Core  Elements  are  used) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-10 

EE  ( I ) 

Elastic  Modulus 

(1) 

11-20 

PR(I) 

Poisson's  Ratio 

(1) 

21-30 

DNS  (I) 

Mass  Density 

(2) 

31-40 

YLD(I) 

Equivalent  Stress  at  First 

Yield 

(3) 

41-50 

ALPHA (I) 

Coefficient  of  Thermal 

Expansion 

NOTES : 

(1)  Repeat  Card  1  for  each  isotropic  core  material  to 

be  defined.  The  first  card  defines  material  property 
set  number  1,  the  second  card,  set  number  2,  and  so 
on.  For  low-modulus  core  materials  in  which  the 
transverse  shear  moduli  are  of  primary  importance, 
setting  EE (I)  =  2G  and  PR (I)  =  0  is  usually  adequate. 

(2)  Mass  densities  are  entered  in  Force-Length-Time 
units  (weight  density /gravity) . 

20 

(3)  The  default  value  of  YLD(I)  is  1.0  x  10  (elastic 
core) . 
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INPUT  SET  7 

ORTHOTROPIC  SANDWICH  CORE  PROPERTIES 

(Skip  this  Set  if  no  Sandwich  Core  Elements  are  used) 


COL 

DATA 

DESCRIPTION 

1 

MTYPE 

Literal  "A"  -  Flag  for 

Orthotropic  Ma-  '.als  Data 

2-10 

El  (I) 

Elastic  Modulus  in  Direction  1 

11-20 

E2  (I) 

Elastic  Modulus  in  Direction  2 

21-30 

E3  (I) 

Elastic  Modulus  in  Direction  3 

31-40 

G12  (I) 

Shear  Modulus  in  Plane  1-2 

41-50 

G13  (I) 

Shear  Modulus  in  Plane  1-3 

51-60 

G2  3 ( I ) 

Shear  Modulus  in  Plane  2-3 

1-10 

PR12  (I) 

Poisson's  Ratio  in  Plane  1-2 

11-20 

PRl 3 (I) 

Poisson's  Ratio  in  Plane  1-3 

21-30 

PR2  3(1) 

Poisson's  Ratio  in  Plane  2-3 

31-40 

DNS  (I) 

Mass  Density 

41-50 

ALPHA1 (I) 

Coefficient  of  Thermal 

Expansion  in  Direction  1 

51-60 

ALPHA2 (I ) 

Coefficient  of  Thermal 

Expansion  in  Direction  2 

61-70 

ALPHA3 (I) 

Coefficient  of  Thermal 

Expansion  in  Direction  3 

NOTES 


NOTES : 


(1)  Repeat  Cards  1  and  2  for  each  orthotropic  core 

material  to  be  defined.  Note  that  all  orthotropic 
materials  must  be  elastic,  and  are  defined  with 
respect  to  the  principal  directions  of  the  material 


INPUT  SET  8 


ORTHOTROPIC  AXIS  DEFINITIONS  FOR  SANDWICH  CORE  ELEMENTS 

(Skip  this  Set  if  no  Sandwich  Core  Elements  are  used,  or  if 
NAXIS  =  0) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

NODE1 (I) 

Node  Numbers  Defining  Origin 
of  Coordinates 

(1) 

6-10 

NODE2 (I) 

Node  Number  Defining  Material 
Direction  1 

11-15 

NODE 3 (I) 

Node  Number  Defining  one 

Additional  Point  in  the  1-2 
Plane  of  the  Material 

NOTES : 


(1)  Repeat  Card  1  to  define  NAXIS  sets  of  axis  directions. 
Orthotropic  material  axes  are  defined  in  terms  of 
existing  nodes  of  the  finite  element  model  as  shown 
in  Figure  A.l. 


INPUT  SET  9 


ELEMENT  CONNECTIVITY,  SANDWICH  CORE  ELEMENTS 

(Skip  this  Set  if  no  Sandwich  Core  Elements  are  used) 


CARD 

COL 

_ 

DATA 

DESCRIPTION 

. 

NOTES 

1 

1-5 

IEL 

Element  Number 

(1) 

6-10 

IPR 

Material  Property  Set  Number 
for  this  Element  (from 

Input  Set  6  or  7) 

11-15 

IAX 

Orthotropic  Axis  Set  Number 
(from  Input  Set  8) 

- 

16-20 

KGEN 

Node  Increment  for  Element 
Generation 

(2) 

21-25 

N  (1 ) 

Local  Node  Number  1 

(3) 

26-30 

N  (2) 

Local  Node  Number  2 

•  • 

56-60 

N  (8 ) 

•  • 

Local  Node  Number  8 

NOTES : 

(1)  Repeat  Card  1  as  required  to  define  all  Sandwich  Core 
Elements.  Elements  must  be  entered  in  ascending 
order,  for  I EL  =  1, 2, . . . ,NELEM.  A  single  blank  card 
is  used  to  terminate  this  section  of  input. 

(2)  A  nonzero  value  of  KGEN  on  the  second  card  of  a  pair 
causes  intermediate  elements  to  be  generated,  by 
incrementing  N(I)  by  KGEN  for  each  succeeding  element 
More  than  one  element  must  be  generated  to  use  this 
feature . 

(3)  Node  point  ordering  for  the  Sandwich  Core  Element  is 
shown  in  Figure  A. 2. 


INPUT  SET  10 


SIZING  PARAMETERS  FOR  STIFFENER  ELEMENTS 

(Skip  this  Set  if  no  Stiffener  Elements  are  to  be  defined) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

IECODE 

Element  Type  Code?  Enter  the 
Number  "3" 

_ 

6-10 

NMAT 

Number  of  Material  Property 

Sets  (Maximum  of  20) 

- 

NELEM 

Number  of  Elements  of  this 
Element  Type 

- 
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INPUT  SET  II 


STIFFENER  ELEMENT  PROPERTIES 

(Skip  this  Set  if  no  Stiffener  Elements  are  used) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-10 

EE  (I) 

Elastic  Modulus 

(1) 

11-20 

PR  (I) 

Poisson's  Ratio 

(1) 

21-30 

DNS  ( I ) 

Mass  Density 

(2) 

31-40 

YLD ( I ) 

Equivalent  Stress  at  First 

Yield 

(3) 

41-50 

ALPHA (I) 

Coefficient  of  Thermal 

Expansion 

NOTES : 

(1)  Repeat  Card  1  for  each  stiffener  material  to  be 
defined.  The  first  line  entered  in  this  Input  Set 
defines  stiffener  property  set  number  1,  the  second 
line  defines  property  set  2,  and  so  on.  Note  that, 
for  face  sheet  stiffeners  (i.e.,  slender  beams) , 

it  is  appropriate  to  set  the  Poisson's  Ratio,  PR(I), 
to  zero. 

(2)  Mass  densities  are  entered  in  Force-Length-Time 
units  (weight/density/gravity) . 

(3)  The  default  value  of  YLD(I)  is  1.0  x  10^° 

(elastic  material) . 
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INPUT  SET  12 


CONNECTIVITY  FOR  STIFFENER  ELEMENTS 

(Skip  this  Set  if  no  Stiffener  Elements  are  used) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

IEL 

Element  Number 

(1) 

6-10 

(Blank) 

11-15 

IPR 

Material  Property  Set  for  this 
Element 

- 

16-20 

(Blank) 

21-25 

KGEN 

Node  Increment  for  Element 
Generation 

(2) 

26-30 

N  (1 ) 

Local  Node  Number  1 

(3) 

31-35 

N  (2 ) 

Local  Node  Number  2 

36-40 

N  (3) 

Local  Node  Number  3 

41-45 

N  (4 ) 

Local  Node  Number  4 

46-55 

THICK 

Stiffener  Width  or  Depth 

(4) 

NOTES : 

(1)  Repeat  Card  1  as  required  to  define  all  Stiffener  Elements 
Face  sheet  and  full-depth  stiffeners  may  be  intermixed 
freely.  All  elements  must  be  entered  in  ascending 
order,  with  a  single  blank  line  entered  to  terminate 
input. 

(2)  A  nonzero  value  of  KGEN  on  the  second  card  of  a  pair 
causes  intermediate  elements  to  be  generated,  by 
incrementing  N(I)  by  KGEN  for  each  succeeding  element. 

More  than  one  element  must  be  generated  to  use  this 
feature. 

(3)  Node  point  ordering  for  all  Stiffener  Elements  is 
shown  in  Figure  A. 3.  N(l)  through  N(4)  are  the  global 
node  numbers  corresponding  to  positions  1  through  4  in 
the  Figure. 
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(4)  The  Stiffener  Element  depth  is  defined  by  the 

positions  of  its  connected  nodes.  THICK  is  the 
element  dimension  perpendicular  to  the  warped 
plane  defined  by  nodes  N(l)  through  N(4). 
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INPUT  SET  13 


SIZING  PARAMETERS  FOR  THIN  SHELL  ELEMENTS 

(Skip  this  Set  if  no  Thin  Shell  Elements  are  to  be  defined) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

IECODE 

Element  Type  Code;  Enter  the 

Number  "5" 

- 

6-10 

NMAT 

Number  of  Material  Property 

Sets  (Maximum  of  20) 

— 

11-15 

NELEM 

Number  of  Elements  of  this 

Element  Type 
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INPUT  SET  14 


THIN  SHELL  MATERIAL  PROPERTIES 

(Skip  this  Set  if  no  Thin  Shell  Elements  are  used) 


CARD 

- — 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-10 

EE  (I) 

Elastic  Modulus 

(1) 

11-20 

PR(I) 

Poisson's  Ratio 

- 

21-30 

DNS (I) 

Mass  Density 

(2) 

31-40 

ALPHA (I) 

Coefficient  of  Thermal 

Expansion 

— 

NOTES : 


(1)  Repeat  Card  1  for  each  Thin  Shell  material  to  be 
defined. 

(2)  Mass  densities  are  entered  in  Force -Length -Time 
units  (weight  density/gravity) . 
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INPUT  SET  15 


THIN  SHELL  ELEMENT  CONNECTIVITY 

(Skip  this  Set  of  no  Thin  Shell  Elements  are  used) 


CARD 

COL 

DATA 

r— - —  1  - -  1 

DESCRIPTION 

NOTES 

1 

1-5 

IEL 

Element  Number 

(D 

6-10 

IPR 

Material  Property  Set  for  this 
Element 

- 

11-15 

KGEN 

Node  Increment  for  Element 
Generation 

(2) 

16-20 

N  (1 ) 

Local  Node  Number  1 

(3) 

21-25 

• 

N  (2 ) 

Local  Node  Number  2 

•  • 

• 

• 

51-55 

N  (8 ) 

•  • 

•  • 

Local  Node  Number  8 

. 

NOTES : 

(1)  Repeat  Card  1  as  required  to  define  all  Thin  Shell 
Elements.  Elements  must  be  entered  in  ascending 
order  for  IEL  =  1,2,...,  NELEM.  A  single  blank  card 
terminates  this  section  of  input. 

(2)  A  nonzero  value  of  KGEN  on  the  second  card  of  a  pair 
causes  intermediate  elements  to  be  generated,  by 
incrementing  N(I)  by  KGEN  for  each  succeeding  element. 
More  than  one  element  must  be  generated  to  use  this 
feature . 

(3)  Node  point  ordering  for  the  Thin  Shell  Element  is 
shown  in  Figure  A. 4. 
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INPUT  SET  16 


SIZING  PARAMETERS  FOR  3-D  SOLID  ELEMENTS 

(Skip  this  Set  if  no  Solid  Elements  are  to  be  defined) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

IECODE 

Element  Type  Code;  Enter  the 
Number  "7" 

6-10 

NMAT 

Number  of  Material  Property 

Sets 

(1) 

11-15 

NELEM 

Number  of  Elements  of  this 
Element  Type 

- 

16-20 

NAXIS 

Number  of  Orthotropic  Axis 
Definitions 

(2) 

NOTES : 

(1)  NMAT  defines  the  total  number  of  property  sets 
to  be  defined  in  Input  Sets  17  and  18. 

(2)  NAXIS  determines  the  number  of  axis  sets  to 
be  defined  in  Input  Set  19. 


iwwvr  set  17 


tsomoric  properties  for  j-d  solid  elmmts 
(Skip  this  Set  if  no  Sol  14  Bleawnts  trt  uood) 


card 

COL 

OATS 

DUCKirTION 

E9 

1 

1-10 

um 

Elastic  Modulus 

(i) 

11-20 

PR(I) 

Poisson's  Ratio 

- 

21-30 

DUS  ( I ) 

Mass  Density 

(2) 

31-40 

YLD(I) 

Equivalent  Stress  at  First 

Yiald 

(3) 

41-50 

ALPHA ( I ) 

Coefficient  of  Thermal 
Expansion 

- 

(1)  Repeat  Cord  1  for  ooch  inotropic  notarial  to  bo 
dof inod. 

(2)  Noon  don sit ion  nro  onto rod  in  Force-Langth-Tiae  units 
(weight  donsity/grsvity) . 

(3)  Tho  dof suit  value  of  YLO(I)  is  1.0  x  1020 
(elastic  oats rial). 


INPUT  SET  It 

OSTNOTNOPIC  MOP  BET  IBS  POE  1-0  SOLID  EL— NTS 
(Skip  this  Sst  if  no  Solid  Elasanto  srs  usad ) 


COL 

r"  - - 

DATA 

— 

DESCRIPTION 

1 

NTYPE 

Litoral  "A"  -  Flap  for 

orthotropie  Notarial o  Data 

2-10 

El  (I) 

Clastic  Nodulus  in  Diraction  1 

11-20 

*2(1) 

Elastic  Nodulus  in  Diraction  2 

21-30 

E3(I) 

Clastic  Nodulus  in  Diraction  3 

31-40 

G12  (I) 

Shaar  Nodulus  in  Plana  1-2 

41-50 

Gl 3 ( I ) 

Shasr  Nodulus  in  Plana  1-3 

51-60 

G23  (I ) 

Shaar  Nodulus  in  Plana  2-3 

1-10 

PR12 (I) 

Poisson's  Ratio  in  Plana  1-2 

11-20 

PR1 3(1) 

Poisson's  Ratio  in  Plana  1-3 

21-30 

PR2  3  (I ) 

Poisson's  Ratio  in  Plana  2-3 

31-40 

DNS  ( I ) 

Mass  Dansity 

41-50 

ALPHA1 ( I ) 

Coefficiant  of  Thansal  Expansion 
in  Diraction  1 

51-60 

ALPHA2 ( I ) 

Coafficiant  of  Thansal  Expansion 
in  Diraction  2 

61-70 

ALPHA 3 ( I ) 

Coafficiant  of  Thansal  Expansion 
in  Diraction  3 

NOTES: 


(1)  Rapaat  Cards  1  and  2  for  sach  orthotropic  solid 

aatsrial  to  ba  dafinad.  Not a  that  all  orthotropic 
Notarial*  must  ba  alastic,  and  ara  dafinad  with 
raspsct  to  tha  principal  diractions  of  tha  Notarial. 


rial. 


IHWT  ir  19 


OKTWOmOtlC  AXIS  DEFINITIONS  ro*  3-D  SOLID  BLMBMTt 

(Skip  this  Set  if  no  Solid  BlssMnts  art  used ,  or  if  MAXIS  -  0) 


COL 

DATA 

DCSCKIPTION 

ES9 

1 

1-5 

NOOBKX) 

Nods  Number  Defining  Origin 
of  Coordinates 

a) 

* 

i 

►- 

o 

NOOB2 (I) 

Node  Number  Defining  Material 
Direction  1 

11-15 

NODE 3 (I) 

Node  Number  Defining  one 
Additional  Point  in  the  1-2 
Plane  of  the  Material 

NOTES t 


(1)  Rspost  Card  1  to  dofins  MAXIS  ssts  of  axis  directions. 
Orthotropic  material  axes  ara  defined  in  terms  of 
existing  nodes  of  the  finite  element  model  as  shown 
in  Figure  A.l. 
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input  trr  20 


COMMSCTXVITY  PON  3-D  SOLID  BLSMNTS 

(Skip  this  Sst  if  no  Solid  BlssMnts  srs  ussd) 


9 


f 


MORS  i 

(1)  MpNt  Cards  1  and  2  aa  raquirad  to  dafina  all  Solid 
11— ant  a.  Bla—nta  nuat  ba  antarad  in  aseandi  ng 
ordar  for  IBL  -  1,2,...,  KBUN.  A  aingla  blank 
card  ia  uaad  to  tarninata  this  saetion  of  input. 

(2)  Aooaptabla  valuaa  for  JWT  ara  I  NT  ■  2  (2  by  2  by  2 
Gauaaian  intagration)  and  I IT  >3  (3  by  3  by  3 
Oauasian  intagration) . 

(3)  A  nonsaro  valua  of  KGEM  on  tha  aacond  card  of  a 
pair  eauaaa  intarmadiata  alaawnta  to  ba  ganaratad, 
by  incrananting  H(I)  by  KQBM  for  aach  auccaading 
al— ant.  Mora  than  ona  alanant  aniat  ba  ganaratad 
to  uaa  thia  faatura. 

(4)  No do  point  ordoring  for  tha  Solid  Blamant  ia  ahoim 
in  Pigura  A. 5.  Noto  that  nodaa  1-8  ara  alwaya 
raquirad,  but  nodaa  9-20  — y  aach  ba  includad  or 
oadttad  aa  daairad. 
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INPUT  SIT  21 

BLmnrr  input  terminator 
(Required  for  all  analysaa) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-80 

(blank) 

(1) 

NOTES t 

(1)  This  Input  lina  terminates  element  input  for  the 
problem  being  defined,  and  ia  required  for  all 
analysis  options. 


INPUT  SET  22 


SIZING  PARAMETERS  FOR  BOUNDARY  CONDITION  DATA 
(Required  for  all  analyses) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

NB  (1) 

Number  of  Type  1  Boundary 

Conditions 

4 

(1) 

6-10 

NB  (2 ) 

Number  of  Type  2  Boundary 

Conditions 

(2) 

NOTES: 


(1)  Type  1  and  Type  2  boundary  conditions  are  defined 
in  Input  Sets  23  and  24,  respectively. 


INPUT  SET  23 


BOUNDARY  CONDITIONS,  TYPE  1 

(Skip  this  Set  if  NB(1)  ■  0  in  Input  Set  22) 


Boundary  conditions  entered  in  this  Input  Set  are 
used  to  constrain  a  range  of  nodes  in  specified  coordinate 
directions. 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

Nl 

Beginning  Node  Number 

(1) 

6-10 

N2 

Ending  Node  Number 

- 

11-15 

INCR 

Node  Number  Increment 

- 

2 

1-5 

JD(1) 

First  Direction  Constrained 

(1) 

6-10 

JD  (2 ) 

Second  Direction  Constrained 

- 

11-15 

JD(3) 

Third  Direction  Constrained 

- 

NOTES : 

(1)  Repeat  Cards  1  and  2  NB(1)  times  to  define  all 

Type  1  boundary  conditions.  Nodes  Nl,  N1  +  INCR, 
Nl  +  2*INCR,...  '72  will  be  constrained  in  the 
directions  JD(I);I  ■  1,2,3.  One  or  more  of  the 
JD(I)  may  be  zero.  Valid  nonzero  entries  for 
JD(I)  are  1,2,3  for  constraints  in  the  x,y  and 
z  directions  respectively. 
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INPUT  SET  24 


BOUNDARY  CONDITIONS,  TYPE  2 

(Skip  this  Set  if  NB(2)  *  0  in  Input  Set  22) 


Boundary  conditions  entered  in  this  Input  Set  are  used 
to  constrain  an  arbitrary  list  of  nodes  in  specified  coordinate 
directions. 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

JD{1) 

First  Direction  Constrained 

(1) 

6-10 

JD(2) 

Second  Direction  Constrained 

- 

11-15 

JD(3) 

Third  Direction  Constrained 

- 

2 

1-5 

ND  (1 ) 

First  Node  Constrained 

(1) 

6-10 

• 

ND  (2 ) 

• 

Second  Node  Constrained 

•  • 

• 

• 

46-50 

• 

• 

ND  (10) 

•  • 

•  • 

Tenth  Node  Constrained 

NOTES : 

(1)  Repeat  Cards  1  and  2  NB(2)  times  to  define  all 
Type  2  boundary  conditions.  The  nodes  ND(I); 

I  =  1,...,10  will  each  be  constrained  in  the 
directions  specified  by  the  JD(I).  Some  of  the 
directions  (JD's)  and  some  of  the  nodes  (ND’s) 
may  be  zero. 
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INPUT  SET  25 


NONLINEAR  LOADING  FUNCTIONS 
(Required  for  nonlinear  analyses  only) 


Loading  and/or  temperature  variations  are  defined  in  this 
Set  as  functions  of  a  single  independent  loading  parameter,  which 
is  continually  increasing  during  the  solution. 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

NFUNCT 

Number  of  Loading  Functions 
to  be  Defined 

(1) 

2 

1-5 

IFUNCT 

Load  Function  Identification 
Number 

(2) 

6-10 

NPTS 

Number  of  Point  Pairs  [X(I), 

Y(I) ]  Used  to  Define  the 

Curve 

(3) 

11-80 

TITL 

Optional  Alphanumeric  Title 

3 

1-10 

X  (1) 

Abscissa  for  First  Data  Point 

- 

11-20 

• 

Y(l) 

• 

Ordinate  for  First  Data  Point 

•  •  • 

. 

61-70 

• 

X  (4 ) 

•  •  • 

•  •  • 

Abscissa  for  Fourth  Data  Point 

71-80 

Y  (4 ) 

Ordinate  for  Fourth  Data  Point 

n 

1-80 

Additional  Pairs  (X(I),Y(I))  as 
shown  for  Card  3 

(4) 
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NOTES: 


(1)  More  than  one  loading  function  will  be  needed  only 
if  non-proportional  loads  are  to  be  considered  in 
the  analysis.  A  maximum  of  fifty  loading  functions 
may  be  defined. 

(2)  Loading  functions  must  be  numbered  between  1  and  50 
inclusive.  The  functions  need  not  be  numbered 
sequentially,  nor  must  they  be  entered  in  any 
particular  order. 

(3)  Up  to  20  data  points  may  be  used  to  define  each  load 
function.  Each  data  point  is  defined  by  an  abscissa 
(the  loading  parameter  value  X(I))  and  an  ordinate 
(the  function  value  Y(I)). 

(4)  Cards  4-7  are  required  for  a  given  loading  function 
only  if  the  number  of  points  exceeds  four.  Unneeded 
cards  should  not  be  entered. 
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INPUT  SET  26 


APPLIED  LOADING 

(Required  for  linear  and  nonlinear  static  analyses) 


CARD 

COL 

DATA 

DESCRIPTION 

NOTES 

1 

1-5 

NODE 

Node  at  which  Load  is  Applied 

(1) 

6-10 

IDIR 

Direction  of  Loading 
(1-X,  2-Y,  3-Z) 

- 

11-15 

NCURV 

Identification  Number  of 

Loading  Function  Describing 
the  Variation  of  the  Load 

(2) 

16-25 

SCALE 

Scale  Factor  or  Load  Magnitude 

(3) 

NOTES : 

(1)  Repeat  Card  1  as  required  to  define  all  applied 
forces  acting  on  the  model.  Loads  may  be  entered 
in  any  order.  Input  in  this  section  is  terminated 
by  a  single  blank  card. 

(2)  The  history  of  each  load  during  the  solution  is 
determined  by  one  of  the  loading  functions  entered 
in  Input  Set  25.  The  parameter  NCURV  is  ignored 
in  linear  analysis. 

(3)  In  linear  analysis,  SCALE  is  the  actual  magnitude 
of  force  to  be  applied  at  the  given  node  in  the 
specified  direction.  For  nonlinear  problems, 
the  value  obtained  from  the  specified  loading 
function  (NCURV)  at  any  stage  of  the  solution  will 
be  multiplied  by  the  value  of  SCALE. 
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